
function [ product_seq_final ] = interchange( products_init_seq, part_names, part_smoothing_intervals, part_weights, products_input_seq, input_part_names, input_smoothing_intervals, input_part_weights )
%interchange Summary of this function goes here
%   Detailed explanation goes here
    %% Initialization
    v = compute_v(products_init_seq, part_smoothing_intervals, part_weights);
    
    p = 1;
    m = size(products_init_seq, 1);
    n = length(part_weights);

    % init loop temporaries
    [curr_v, curr_part_names, curr_smooth_intervals, curr_part_wts] = sort_v(v, part_names, part_smoothing_intervals, part_weights);
    curr_prod_seq = products_init_seq;

    max_passes = 1000;
    for k = 1:max_passes
        interchange_count = 0;
        %% propose interchange
        for i = 1:m-1
            for j = i+1:m
                % does interchange [i,j] reduce v(p)?
                candidate_seq = swap_products(i, j, curr_prod_seq);
                if not(is_sequence_feasible(candidate_seq))
                    continue
                end

                candidate_v = compute_v(candidate_seq, curr_smooth_intervals, curr_part_wts);
                if candidate_v(p) < curr_v(p)
                    %% evaluate interchange
                    % will this interchange increase v(h) for at least one h < p
                    if p > 1 && sum(candidate_v(1:p-1) > curr_v(1:p-1)) > 0
                        continue; % reject interchange
                    end

                    % will this interchange cause some v(h) become larger tha v(p) for any h > p
                    if p < n && sum(candidate_v(p+1:n) >= curr_v(p)) > 0
                        continue; % reject interchange
                    end

                    i
                    j
                    curr_prod_seq = candidate_seq; % execute interchange
                    interchange_count = interchange_count + 1;
                    
                    % update the list
                    [curr_v, curr_part_names, curr_smooth_intervals, curr_part_wts] = sort_v(candidate_v, curr_part_names, curr_smooth_intervals, curr_part_wts);
                end
            end
        end
        
        if p < n
            p
            interchange_count
            p = p + 1;
        else
            % termination condition
            if interchange_count == 0
                curr_v
                curr_part_names
                k
                product_seq_final = candidate_seq;
                disp_result(products_input_seq, input_part_names, input_smoothing_intervals, input_part_weights, product_seq_final, curr_part_names, curr_smooth_intervals, curr_part_wts);
                break;
            end
            interchange_count
            p = 1;
        end
    end
end

function [] = disp_result(seq_input, input_part_names, input_smoothing_intervals, input_part_weights, seq_final, final_part_names, final_smoothing_intervals, final_part_weights)
% plot final q_ik's vs preprocessed q_ik's

    [~, q_input] = compute_v(seq_input, input_smoothing_intervals, input_part_weights);
    [~, q_final] = compute_v(seq_final, final_smoothing_intervals, final_part_weights);
    
    n = length(input_part_names);
    figure % create new figure
    m = size(seq_input, 1);
    x = linspace(1, m, m);
    for i = 1:n
        part = input_part_names(i)

        subplot(3,2,i); % first subplot
        find(final_part_names == part)
        plot(x, q_input(part,:), '--ro', x, q_final(find(final_part_names == part),:), '--go', x, ones(m, 1)*mean(q_input(part,:), 2), '-r', x, ones(m, 1)*mean(q_final(find(final_part_names == part),:), 2), '-g');
        title(sprintf('Part %d', part));
    end
    %plot()
end

function [ is_feasible ] = is_sequence_feasible(sequence)
% is_sequence_feasible: check if the product sequence satisfies
% earliest release and due date constraints
    is_feasible = true;
    m = size(sequence, 1); % number of products
    
    for i = 1:m
        if i < cell2mat(sequence(i,2)) || i > cell2mat(sequence(i,3))
            is_feasible = false;
            break;
        end
    end
end

function [ v_sorted, part_names_sorted, part_smoothing_interval_srtd, part_wts_srtd ] = sort_v(v, part_names, part_smoothing_intervals, part_weights)
%% sort list of deviations and do other housekeeping
    % sort v's
    [v_sorted, sort_index] = sort(v, 'descend');
    % sort part-related data to be consistent
    part_names_sorted = part_names(sort_index);
    part_smoothing_interval_srtd = part_smoothing_intervals(sort_index);
    part_wts_srtd = part_weights(sort_index);
end

function [ prod_seq ] = swap_products(i, j, tmp_prod_seq)
%% return a copy of tmp_prod_seq with rows i and j swapped
    prod_seq = tmp_prod_seq;
    prod_seq(i,:) = tmp_prod_seq(j,:);
    prod_seq(j,:) = tmp_prod_seq(i,:);
end

%% compute v's
function [ v, q ] = compute_v( product_seq, smoothing_intervals, weights )
    q = compute_q(product_seq, smoothing_intervals);
    q_i = max(q, [], 2);
    q_i_bar = mean(q, 2);
    v = weights'.*((q_i - q_i_bar)./q_i_bar);
end

%% compute q's
function [ q ] = compute_q( product_seq, smoothing_intervals )

    m = size(product_seq, 1); % number of products
    n = size(cell2mat(product_seq(1,4)), 1); % number of parts
    q = zeros(n,m);
    r = q;
    for i = 1:n
        %part_name = part_names(i)
        r(i,:) = get_part_req(i, product_seq);
        q(i,:) = get_q_ik(r(i,:), smoothing_intervals(i));
    end
end

%% get part requiremnt sequence for part at index part_index
function [r] = get_part_req( part_index, product_seq )

    m = size(product_seq, 1); % number of products
    r = zeros(m,1);
    for i = 1:m
        t = cell2mat(product_seq(i, 4));
        r(i) = t(part_index);
    end
end

function [q_ik] = get_q_ik( r_i, smoothing_interval )
    
    q_ik = r_i;
    m = length(r_i);
    for j = 1:m
        if smoothing_interval > 1
            for k = 2:smoothing_interval
                idx = mod((j + k-1), m);
                if idx == 0
                   idx = m;
                end
                q_ik(j) = q_ik(j) + r_i(idx);
            end
        end
    end
end